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Abstract — This paper examines filtering on a sphere, by first 
examining the roles of spherical harmonic magnitude and phase. 
We show that phase is more important than magnitude in 
determining the structure of a spherical function. We examine 
the properties of linear phase shifts in the spherical harmonic 
domain, which suggest a mechanism for constructing finite- 
impulse-response (FIR) filters. We show that those filters have 
desirable properties, such as being associative, mapping spherical 
functions to spherical functions, allowing directional filtering, 
and being defined by relatively simple equations. We provide 
examples of the filters for both spherical and manifold data. 

Index Terms — spherical harmonics, phase, FIR filtering 

I. Introduction 

While the importance of phase information in the ordi- 
nary Fourier transform on Euclidean domains is well- 
understood, it may be argued that a comparable understanding 
of phase does not exist when applied to the sphere S'^. It 
is important to fill that gap because spherical data arise in 
numerous fields, including measurements of atmospheric pres- 
sure, cosmic microwave background, or surface reflectance. 
Consequently, methods for filtering, spatio- spectral and spatio- 
scale analysis of spherical data using spherical harmonics have 
attracted the attention of the research community iTTllfTSllfTSl . 
However, existing research has not provided a detailed under- 
standing of the importance of phase of spherical harmonics, 
either on its own or in relation to the harmonic magnitude. It 
is the purpose of this paper to explore both the properties and 
the applications of phase information for spherical harmonics. 
Our exploration of phase leads to a new method of filtering 
on the sphere which has desirable properties such as allowing 
for directional filtering, associative construction, and mapping 
spherical functions to spherical functions. 

We begin by recalling important properties of phase from 
the ordinary Fourier transform. Let the Fourier transform of 
/ on the real line R be denoted F = TT {f}. The transform 
may be split into magnitude |F| and phase e^^ as follows: 
F = \F\e^^. If we set the phase to zero and invert the trans- 
form, i.e., by letting G = \F\ and g = J^T~^ {G}, then the 
resulting function g must be a positive-definite function with 
its maximum value at the origin as illustrated in Figure [ij^)- 
In fact, it is well-known that g is the autocorrelation of the 
function obtained from the inverse Fourier transform of the 
square-root magnitude, V"]^. Recall that the autocorrelation 
is the integral 

«/(^) = / r{x)f{x^s)dx. (1) 
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Fig. 1: In part (a), the inverse Fourier transform of the 
magnitude-only spectrum from a white noise sample is shown. 
The result is a positive definite function, which is an autocorre- 
lation. Part (b) shows the analogous result for the sphere, the 
inverse spherical harmonic expansion of the magnitude-only 
spectrum as discussed in this paper. 

Then, it is well known that Af = J^T {a/} = FF* = 
Hence, if its phase is set to zero, each function is turned into 
an auto-correlation. Similarly, if we operate on photographic 
images, and swap the phase of image / with that of image 
g, i.e., if we let / = J^l ^ ||F|y^|, then it is well-known 

1 21 1 that / resembles g in appearance much more than /. 
As we show below, similar demonstrations are possible with 
harmonic analysis on the sphere. Figure [TJb) gives a preview 
by showing the result for S'^ that corresponds to removing 
phase entirely, using methods described later in this paper. 

Given the importance of phase in determining the appear- 
ance and structure of data, we might expect that phase infor- 
mation is crucial for applications such as viewpoint-invariant 
3-D shape recognition in computer vision. Yet, despite that, 
spherical harmonic invariants for shape recognition are in 
fact phase-blind magnitudes |[T4lll23ll . A recent paper ifTSll 
showed that by using phase- sensitive invariants derived from 
the bispectmm of spherical harmonics, both discrimination and 
robustness to noise improve. There has been little comment 
otherwise on the role that phase plays in determining spherical 
functions, which is one of the motivations for this paper. 

There is an extensive literature on filtering on the sphere. 
The mathematical aspects of that literature are discussed later, 
but here we summarize the contributions made to date and 
indicate how our approach is different. A pioneering work by 
Driscoll & Healy | 8 | describes convolution of two functions, 
/ and h, on the sphere by integrating the rotated version 
of h with /. The result is a function g on 5*0 (3), which 
may be parameterized by the Euler angles denoted a, /3, 7 
as ^(a,/3,7). For a fixed 7, the output g^{a^P) is also a 
function on 5^. However, for the fixed 7, there is no harmonic 
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space product expression that is analogous to the famiUar 
product F • G on R. A product expression occurs, however, 
when g is projected back to the sphere by integrating out one 
angle of rotation, which is equivalent to using rotationally- 
symmetric version of the filter. Those filtering methods are 
adopted in various papers |4||7|| 18|| 19||28||29||311|30||33l. 
In this paper, we propose an alternative filtering method 
which has two main differences to the previous methods: it 
uses a relatively simple harmonic-space product expression 
which works for directional filters, and it provides a means 
for constructing finite-impulse response (FIR) filters. There 
are two key insights which lead to our results: the use of 
vector-matrix expressions that simplify construction of filters 
in the harmonic domain, and the analysis of phase in spherical 
harmonics. 



II. Background material 

The theory of 3-D rotations, group representations, and 
spherical harmonics is well-described in various books 
ll5lll9ll[TQl . We review key concepts and results briefly in this 
section. Every rotation is represented by three Euler angles, 
denoted a, f3, and 7, with a e [0, 27r] denoting rotation 
about the z axis, (3 G [0,7r] denoting rotation about the y 
axis, and 7 G [0, 27r] also denoting rotation about the z axis. 
Two of the angles, a and /3, parameterize the sphere S'^, 
with < a < 27r representing longitude (angle measured 
counter-clockwise from the positive x-axis in the x-y plane), 
< /3 < TT the colatitude (angle with respect to positive 
2:-axis). Each unit vector u lying on S'^ may be described 
parametrically as 



u = [cos(a) sin(/3), sin(Qf) sin(/3), cos(/3)]^ . 



(2) 



The Laplace spherical harmonics form an orthonormal basis 
for functions on S'^. For each non-negative integer i, and 
integer —£< m < i they have the functional form 

Yp{u) = Yl^{p,a) = Q^Pr(cos/3)e-^'^", -i < m < i. 

(3) 

Here Pp are the associated Legendre functions, which are 
real-valued, and the normalization constants are 



/2f + 1 (t-m)\ 
47r (i^m)\' 



(4) 



For convenience, we denote the integral on the sphere as 
follows 



f{u)du ■■ 



f{P,a)sm{P)dpda. (5) 



With the above definitions, we have orthonormality 

Yr{u)Yp^{uydu = 5,p5mn- 
Defining the £, m-th spherical harmonic coefficient as 



f{u)Yr{urdu, 



(6) 



(7) 



every square-integrable function / on the sphere may be 
expanded as follows: 



f{u) = E 



F^Yl 



(8) 



i=0 m=-i 



It proves very convenient to write ([8]) as a series of vector 
inner products. The vector notation not only reduces the clutter 
of indices but also provides linear- algebraic insights which 
led to the main results in this paper. First, let us combine all 
coefficients for a given i into a 1 x (2^ + 1) row vector 



(9) 



We refer to F£ below as the ^-th Fourier coefficient of /. 
Second, let Y^ denote the (2^+1) x 1 column vector containing 
the spherical harmonics at the ^-th frequency 

Ye{u) = [Yf'iu), . . . , . . . , Yf{u)] ^ . (10) 

Using these vectors, we write ^ as 



(11) 



i=0 



This relatively- simple equation is, as we see below, the key 
to understanding the role of phase and the implementation of 
filtering. Suppressing indices, we may "visualize" the equation 
as a sum of inner products of progressively longer odd-length 
vectors: 



f{u) = FF+f FFF ] 



Y 
Y 
Y 



-[ FFFFF ] 



Y 
Y 
Y 
Y 
Y 



(12) 

Here the symbols F and Y inside the vectors denote appro- 
priate elements of the F^ and Y^ vectors with their indices 
suppressed. The same expansion ^TT) , using only a finite 
bandwidth L is 



L-l 



/^(ii) = EF,r,(^). 



(13) 



i=0 



A. Magnitude and phase 

It is interesting at this point to consider what might be 
the "phase" of spherical harmonic coefficients. We may, for 
example, split each scalar coefficient Fp into magnitude 
\Fp\ and "phase" e^^^^'^) = Fp/\Fp\. However, that phase 
does not behave analogously to the phase spectrum of the 
ordinary Fourier transform. When a function / on R is 
translated by x x-\-t, then its phase spectrum transforms as 
(/) + 27nyt, where u is the frequency. There is no similar 
behavior for the scalar "phase" (j){£^m) of spherical harmonic 
coefficients just defined. 

To determine a more appropriate interpretation of phase for 
the spherical harmonics, we examine their rotation property. 
Every 3-D rotation is represented by 3 x 3 orthogonal matrix R 
with determinant +1. If, for some R, we have g{u) = f{Ru) 
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for all u e S'^, so that ^ is a rotated version of /, then the 
Fourier coefficients of g are related to those of / as follows: 



G, = F,D,{R). 



(14) 



Here Di{R) is a {21 + 1) -dimensional unitary matrix, known 
as the Wigner D-matrix, with the homomorphic property 
Di{RS) = Di{R)Di{S) for every pair of rotations R, S. The 
elements of the /^-matrices are separable in the Euler angles: 

L>P(a,/3,7) = e-^'^"d7^^(/3)e-^'^^, -i<m,n<i. 

(15) 

Here d'^'^ is the "little" Wigner d-f unction, which is real- 
valued for the z-y-z Euler angles as defined earlier. In par- 
ticular, the spherical harmonics are, up to a scale factor, the 
elements of the middle column of the corresponding D matrix: 

(/^, o^) = ^DY'^a, /3, 0), < m < (16) 
where, from (|4]), we have the constants 



2^- 



47r 



(17) 



We now obtain an appropriate definition of phase. Each co- 
efficient vector may be separated into magnitude \\Fi\\ and 
"phase" determined by the unit-length vector = F^/HF^H, 
so that Fi = 1 1 1 1 [7^. Under a rotation of the underlying 
function on the sphere, each coefficient F^ transforms by 
([14]). Since the Wigner matrices are unitary, we see that 
the magnitude \\Fi\\ will not change, but the phase vector 
transforms as /J^ UiDi{R). That behavior matches the 
linear shift property of the ordinary Fourier transform under 
translation. 

in. Properties of magnitude and phase 

Any attempt to define magnitude and phase is made clearer 
with an exploration of properties. We proceed by looking 
to the familiar properties on the real line R. There, the 
Fourier transform of the Dirac delta function is constant for 
all frequencies: 



/CO 
5{x)e-^'^'-'dx = 1. 
-OO 



(18) 



When the delta function moves to any location xq, only the 
phase of F changes: F{uj) = e~^^^^. 

Unlike R, it is not obvious how to define a delta function 
on S'^, since the sphere does not have an "origin" with unique 
properties. We take our inspiration from the equation 

/OO 1 
^-jux^juyj^^ (19) 
-OO 27r 

Following [24, pp 594-595], we define for any two points u, 
w on S'^, the function (with f denoting conjugate-transpose) 



£=0 



It is easy to show, for any function /, that 

S{u — w)f{u)du — f{w). 



/ 



(20) 



(21) 
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Fig. 2: The magnitude spectrum of the S function on the sphere 
5*^ . It increases with frequency £ unlike the counterpart on the 
real line R. 



Suppose now that we have a delta function S{u — n) located 
at the North pole n = [0,0, 1]^ of the sphere. We take the 
angular coordinates of the north pole to be a) = (0,0). 
The spherical harmonic coefficients are 



S{u-n)Yp{uydu = Yl^{0,Oy =cfP^{l). (22) 



Since P7^(l) = 1 for m = 0, and zero for other m, we have 
that the Fourier coefficients of the delta function are 



if m = 
otherwise 



(23) 



Consequently, the vectors F^ are all zero except in the middle 
entry, which is equal to q. The magnitude spectrum of the 
delta function is independent of its position on the sphere since 
it is rotation invariant. Therefore, for a delta function located 
anywhere on the sphere, we have 



\FA 



2i- 



(24) 



This result has been published before |[22| eq 64]. Therefore, 
unlike on the real-line, the magnitude spectrum of a delta 
function on S'^ is not constant, but rather increases with 
frequency £ in the manner illustrated by Figure [2] 

A second example of spherical harmonic expansion is 
obtained for the Fisher Von-Mises distribution ifTTl . which is 
defined with scale parameter k, and mean /i G 5*^ as 



47rsinhA>: 



(25) 



Let 11 = the north pole. Then u = cos /3, and the distri- 
bution may be expanded in terms of the Legendre polynomials 
as follows 1251 eq 20]: 

lf:(2£ + l)^^i/^P,(cos/^) (26) 
47r f-^ /i/2(t) 



f{u;n,K) 
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(a) 




(b) 

Fig. 3: The Fisher- Von Mises distribution is shown in part (a) 
for /i: = (green), k, = 2.b (blue), z^: = 10 (red). The polar 
angle (3 is shown from — tt to tt to center the distribution, 
though in practice only [0,7r] is used. Part (b) shows the 
magnitude spectrum for the corresponding values of hz. 



Here, /j^ is the modified Bessel function of order u. It is easily 
shown from this formula that FI^ = unless m = 0, and that 



(27) 



Figure [Sja) and (b) show the distribution and its magnitude 
distribution for various k,, 

The magnitude-phase partition of spherical harmonic coef- 
ficient vectors, Fi = invites us to consider which 
is more important to the structure - magnitude or phase? 
Considering that F^ is 2^ + 1 dimensional, but \\F£\\ is a 
scalar, we should expect that the magnitude constrains only 
a small portion of the structure: phase should be much more 
important. Suppose that / is a real-valued function. Then, 
due to conjugate symmetry, F^^ = ( — 1)^ (F^^)*, so that 
each Fi vector has only 2^ + 1 real degrees of freedom. Of 



these, one degree of freedom is constrained by the magnitude 
IIF^II, leaving 2£ degrees of freedom in the phase vector 
Ui = Consequently, for a function / defined by 

its first £ < L spherical harmonic vectors, the magnitude 
spectrum HF^H for £ < L constrains only L + 1 of the total 



^2^ + l = (L + l)2 



(28) 



i=0 



degrees of freedom. This means that the phase spectrum, 
{/7o,/7i,...,/7l} constrains the remaining L(L + 1) degrees 
of freedom, or as a percentage 



100- 



(29) 



For example, for L = 10, the phase spectrum constrains 
roughly 91% of the signal. Similar estimates can be derived 
for complex-valued signal. This is unlike the situation on R, 
where the phase and magnitude each constrain exactly half the 
degrees of freedom of the transform coefficients. We see that 
on the sphere, phase is much more important to a signal than 
magnitude. 

We now describe an important property of magnitude, that it 
only determines an axially- symmetric component of a function 
on S'^. First, suppose that / is axially-symmetric about the 
North pole n = [0,0, 1], i.e., f{Ru) = f{u) for all rotations 
R such that Rn = n. Then, from ([14]), we see that its spherical 
harmonic coefficient vectors satisfy F^ = FiDi{R) for such 
rotations R. From ( p3] ) we see that this is only possible if 
F^ = for m ^ 0, i.e., Fi has only one non-zero element, 
which located at the center of the 2£ -\- 1 dimensional row 
vector. Define the 2£ -\- 1 row vector Qi with 1 in its middle 
entry: 

Q^ = [0,...,0,1,0,...,0] (30) 

Therefore, functions axially-symmetric about the North pole 
satisfy Fi = F^Qi. Consider now a spherical function / 
about whose spherical harmonic coefficient vectors {F^}^o 
we know only the magnitude spectrum \\Fi\\. Although there 
are infinitely many functions with the same magnitude spec- 
trum, we see that one possibility is the axially-symmetric 
function whose spherical harmonic coefficient vectors satisfy 
= IIF^IIQ^. Since such functions are entirely determined 
by their magnitude spectra, we conclude that the magnitude 
component does not determine any more than an axially- 
symmetric version of the function. 



A. The view from 6'0(3) 

The properties of spherical harmonics become clearer when 
S'^ is expressed as a homogeneous space for the group 50 (3) 
of 3-D rotations. That group consists of all 3 x 3 orthogonal 
matrices with determinant +1. Every function / whose domain 
is S'^ may be "lifted" to a corresponding function / on 50(3) 
by the mapping f(R)= f{Rn), where n is the North pole as 
before. _ 

Note that / is constant on rotations that fix n, and, therefore. 



/(a,/3,7)=/(a,/3,0), V7. 



(31) 
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The elements {1^7^^}^,m,n of the Wigner I)-matrices men- 
tioned above form an orthogonal basis for functions defined on 
S0{3). We review some basic facts; see |16||27| for details. 
The Fourier transform on 50(3) consists of matrix-valued 
coefficients, one at each "frequency" £, where £ = 0, 1, — 
The ^-th Fourier coefficient is 



F{£) 



50(3) 



f{R)De{R)UR. 



(32) 



Here the symbol f denotes matrix hermitian transpose and 
dR = {87r'^)~^ sm{P)dadf3 is the normalized, rotation- 
invariant, measure on SO (3). A 3-D rotation of /, which 
equivalently translates the function on SO (3), produces a 
corresponding transformation of the Fourier coefficients: 



g{R) = f{SR) ^ G{i) = F{£)D,{S). 



(33) 



With ( |32| ), the function / may be expanded as follows 

CO 

f{R) = ^{2£ + l)Trace [F{£)D,{R)] . (34) 



£=0 



The factor (2^ + 1) corrects for the L2 norm of the matrix 
elements: 

1 



SO{3) 



\D^''{R)\UR: 



2^ + 1' 



(35) 



If / = /, or in other words / is a function lifted from 5'^ 
as described above, then from eqns. ( p3] ) and ( [3T] ) we obtain 
that only the middle row of the {2£ + 1) -dimensional Fourier 
coefficient matrix F{£) is non-zero. 



m^0=> F{£f 



0. 



(36) 



Moreover, from ([16]), we obtain (see 1 13|) that for lifted func- 
tions, the middle row of each of the 5*0(3) Fourier coefficients 
contains, up to a scale factor, the Fourier coefficients of the 
spherical harmonic expansion: 



F^£r = ^Fr. 



(37) 



It is easy to verify that, with ( [37] ) and (16), the expansion ( [34] ) 
reduces to ([5]). 

IV. Filtering 

Filtering, which is convolving a signal / with a kernel 
h, is the most basic signal processing operation. Due to 
the non-commutative nature of the rotation group 50(3), 
filtering is not easy to define on the sphere 5^ in terms of 
an integration operation. In this section, we review several 
plausible definitions of filtering on the sphere, including a 
novel approach using lifting as defined in Section |III-A[ Any 
useful definition of filtering, denoted g = f h, has the 
following properties. 

Fl It maps spherical function / to an output g which is 

also a spherical function; 
F2 The kernel h may be either axially- symmetric or 

directionally selective; 
F3 It is associative, so that {f^hi)i<h2 = /^(/^-i^o ^2), 
where is a suitable convolution on the domain of 
the filters hi and h2\ 



F4 



It permits of FIR filtering, which is defined as a finite 
weighted sum of rotated copies of the function 



N-l 



(38) 



k=0 



The weights {b^} and the rotations {Rk} define the 
operation of the filter. 

Note that simply rotating a function on 5^ is a filtering 
operation, and is a special case with = 1 of FIR filtering 
described in ( [38] ). 

We describe three approaches to filtering, and verify for 
each approach whether it has properties F1-F4. The first 
approach, which we call the "rotation" approach, to filtering 
is to define convolution on 5^ analogously to that on the real 
line, so that the domain is 5^ and the translational variable, 
which in this case is rotation, is applied to the kernel 



g{R) = {f^h){R) 



52 



P{u)h{Ru)du. 



(39) 



This definition produces a function g on 50(3), not 5^, and 



thus does not have property Fl above. Although (39) is not 
suitable, it is worth noting that it has a convenient expression 
in terms in the frequency domain. The Fourier coefficients of 
G on 50(3) are outer-products of F^ and Hi as follows: 



G{£) = ^fIh, 



(40) 



with constants r]i = 1/(2^ + 1). Note that both sides of the 
above are 2£ -\- 1 dimensional matrices. A proof using vector 
notation is provided in Appendix A; a similar result using 
indices is shown by Kostelec & Rockmore |16|. 

The second definition of filtering, which may be called "left 
convolution", is described and analyzed in Driscoll & Healy 
|8|. For /, h in 1/^(5^), and n the North pole, the convolution 
is obtained by integrating over rotations: 

g{u) = {fi.h){u) = [ f{Rn)h{R-^u)dR (41) 

It has been shown |8| that spherical harmonic coefficients of 
g are obtained by 



977- 



(42) 



The coefficient is defined in ( 17 ). Note that only the central 



coefficient of Hi vector is involved in producing the output 
vector. Therefore, the phase of Gi is, up to sign, the phase of 
Fl, so that filtering defined by ( [4T] ) essentially preserves phase. 
Since only the central coefficient is used, this definition only 
works for axially- symmetric kernels, and therefore does not 
have property F2. 

The third definition of filtering, which is our proposal, is to 
convolve both functions / and h on 50(3), and project the 
result back to the sphere. The basis of this definition is the 
convolution of two functions on 50(3), which is defined as 



9{V) = I 
Js 



h{R)f{R-^V)dR. 



(43) 



SO{3) 
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Using ([32]), it may be shown (see Appendix B) that the Fourier 
transform of g is related to those of /, as follows: 



G{1) = F{£)H{£) 



(44) 



We define the projected convolution on 6*0 (3) as follows. Let 
/ be a function on S'^, and h a kernel function on 5*0 (3), 
whose construction is discussed later. Define 



9{u) Ps2 



h{R)f{R-^V)dR 



50(3) 



(45) 



Here the projection operator P52 converts the three variable 
function on SO (3) to a two- variable function on S'^ by 
integration. Since a, /3 are sufficient to parameterize the sphere 
by ([2]), we integrate over the Euler angle 7 as follows: 

Ps^ [9{R)] = ^JJ ^(^' ^)^^- ^46) 

The properties of ( [45] ) are now established. For each i > 0, 
let Gi be the {2£ + 1) -dimensional row vector of spherical 
harmonic coefficients for g, the filter output, Fi is the cor- 
responding vector for /, and H{i) is the square {2i + 1)- 
dimensional Fourier coefficient matrix obtained from ( [32] ) for 
the filter kernel h. 

Theorem 4.1: With the notation as above, the projected 



convolution ( 45 ) has the Fourier representation 



(47) 



Furthermore, it satisfies all four properties F1-F4. 

Proof: Using as defined in (30), we have from ( [37] ) 

that 

m = ^oQjF^ (48) 



Let the integral in ( [45] ) be denoted 



g(y) = [ h{R)f{R-W)dR, 

JSO{3) 



This is a function on SO {3), and its Fourier coefficients, from 
( [44] ) and ([48]) are 



G{i) = ^QjF,H{i) 



(50) 



Consequently, using ( 34 ), and rearranging terms, we see that 
(|45l) becomes 



\ Trace 



dj 
2^ 



Using ( [T5] ) and ( [T6] ), we obtain that 

CO 

g{u) = YF,H{i)Y,{u) 



(51) 



(52) 



i=0 



from which, by comparing with ( [T7] ), the result (47) follows. 

We now verify satisfaction of the properties. Fl is obvi- 
ous. To prove F2, note that the case of axially- symmetric 
filtering, which is provided by ( [42] ), is obtained by setting 
H{£) = ^H^l£, where is the 2£ + 1-dimensional identity 
matrix. Tfie case of directionally-selective kernels is estab- 
lished by the FIR filtering property shown next. The property 



F3 follows from the transform ([47|, since {f i^hi)i<h2 gives 
coefficients F£Hi{£)H2{i), which by (44) is the same result 
for transforming (/ii ^2)- To prove F4, note that the FIR 
filter is easily shown have the Fourier transform 



N-l 



(53) 



k=0 



Examples of directional FIR filters are given below. ■ 
The three types of filtering, including our proposal, are 
summarized, with a "visualization" in terms of vectors and 
matrices for the first coefficient ^ = 1 in Table U 



Equation 


Domain of g 


Illustration of G for = 1 
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SO{3) 




GGG 
GGG 
GGG 


oc 


F 
F 
F 


[ HHH ] 







[ GGG ] oc [ FFF ] H 




45 






[ GGG ] = [ FFF ] 


HHH 
HHH 
HHH 





TABLE I: Summary of filtering methods g = / ^ on the 
sphere, with the proposed method on the bottom row. In the 
rightmost column, G, F, and H represent generic (^indices 
suppressed) harmonic coefficients of output, signal, and kernel, 
respectively. 



Table [I] helps to explain why ( [45] ), our proposed method 
of filtering, is a "phase- sensitive" method. The previously- 
known method of filtering on 5^, (41), allows only axially- 



symmetric filters to be used. Axially-symmetric filters, as 
shown in Section [III] are determined essentially by their 
magnitude spectrum and therefore have no phase component. 
Furthermore, the effect of the previous method R^ on the 



(49) 

input function / only changes its magnitude components (up to 



sign) and not its phase. In contrast, the proposed method (45 ) 
allows freedom in designing the phase response of the filter, 
and consequently provides a wide range of filtering effects 
on both the magnitude and phase components of the input 
function. 

Here are some examples illustrating the projected convolu- 
tion approach. The simplest filter is rotation: g{u) = f{Ru) 
for some rotation R. The coefficients of this filter are uni- 
tary: H{£) = Di{R) for all I. Therefore, the phase of the 
signal coefficients, Ui = Fi/\\F£\\, is transformed linearly: 
Ui ^ UiDi{R). A second, relatively simple filter, is defined 
on S'^ and has an axially-symmetric kernel. As discussed in 
the proof, its S0{3) Fourier coefficients are H{t) = ^H^h, 
where Ii is the 2£ + 1-dimensional identity matrix and are 
its central spherical harmonic coefficients.. 

A third example is a local average of f{u) and two 
neighbors at lines of higher and lower latitude, respectively. 
This 3-point filter may be described as 

g{u) = O.bfiu) + 0.2b f{R^,u) + 0.2bf{Rj^u). (54) 

Here, Rf^^ is the rotation that moves the North pole n to 
the axis [sin/3o, 0, cos /3o], which is aligned along the a = 
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meridian. The corresponding filter has the Fourier transform 

H{i) = 0.5Ie + 0.25De{Ro) + 0.25De{Ry . (55) 

This is an example of a filter that is not axially symmetric, and 
which may be implemented using our projected convolution 
approach. 

A fourth and more detailed example, is obtained by applying 
the associativity property F3. If we apply the FIR filter ( [55] ) 
to the output of an axially symmetric filter (such as the 
Fisher- Von Mises) then the result a weighted combination of 
symmetric filters arranged on the latitude axis. Such a filter 
may be described as 



g = 0.5/ ^ + 0.25/ i^hn^ 0.25 f ^ h^T . 



(56) 



In the frequency domain, the filtering is described as Gi = 
FiHcii), where the cascaded filter has coefficients 

H,{£) = —H^ [0.5/, + 0.25i^,(^/3o) + 0.251^, (i^^^o)^] • 

(57) 

We refer henceforth to the H{£) matrices as the transfer 
function of the filter. An important concept in filtering is the 
impulse response. We obtain an analogous impulse response 
for the projected convolution using the inverse spherical har- 
monic expansion of ( [47] ), where the input coefficient vectors Fi 
are those for the spherical impulse as defined in (23 ). However, 
such an impulse response is incomplete, since it determines 
only a function on S'^ and not the kernel on 6*0 (3). Hence, 
the transfer function is required as the complete description of 
the filter. 

V. Experiments 

We implemented the projected convolution and validated 
its properties in three experiments. In each experiment, we 
use the method of computing spherical harmonics described 
by Chung et al O, and the software accompanying that 
paper. That method is called iterative residual fitting (IRF), 
which progressively computes at each frequency the spher- 
ical harmonic coefficients which best fit, in a least- squares 
sense, the residual obtained by subtracting from the data a 
weighted representation from the next lower frequency l—l. 
Intuitively, the progressive weighted coefficient fitting provides 
a "windowing" effect similar to that known to reduce Gibbs 
phenomenon for Fourier series on the real line. We also 
used data and software accompanying the papers ISl |[T6l 1261 . 
MATLAB code for all of our experiments is available online 
113. 

In our first experiment, we construct a simple 5 -tap FIR 
lowpass filter. Figure [4] shows the reconstruction of an impulse 
defined in the frequency domain using ( 23 ), using bandwidth 
L = 63 as defined in ( [T3] ). The reconstruction is obtained 
using IRF. The filter is defined by the input-output equation 

g{u) = 0.5f{u) + 0A25f{Riu) + 0A25f{Rju) + 

0.125/(i?2^) + 0A25f{Rju) (58) 

The two rotation matrices, Ri and R2 are defined as follows: 
Ri has Euler angles a = 0, (3 = 7r/32, and 7 = 0, while 
R2 has Euler angles a = 7r/2, (3 = 7r/32, 7 = — 7r/2. The 



rotation Ri shifts the North pole down by 7r/32 along the 
0^ (Greenwich) meridian, and the second R2 shifts by the 
same amount along the 90^E meridian. This filter is easily 
implemented in the frequency domain using the coefficient 
matrices 

H{i) = 0.5Ii + 0.125 [D^{Ri) + Di{Ri)^ + 

0.125 [De{R2)^De{R2y] (59) 

Figure [4] shows the impulse response of the filter, illustrating 
four smaller lobes surrounding a central main lobe. The 
impulse response is computed using the inverse transform of 
the projected convolution, which is obtained in the frequency 



domain using ( [47] ). As a simple summary of the frequency 
response, we computed the norm \\H£\\ at each frequency £, 
and plot the results in Figure [4] The plot shows values after 
normalization by the magnitude of the impulse's coefficient 
vectors as shown in ( [24] ). It can be seen that the norm, which 
measures the "size" of the magnitude component, generally 
decreases with £, clearly the characteristic of a lowpass as 
expected from a local average. 

In our second experiment, we filtered a binary world map 
shown in Figure [5] In the same figure, we show a reconstruc- 
tion of the world map using bandwidth L = 63, as in ( 13 ) but 
where the coefficients are obtained by IRF. The value of 
IRF in this case is that it suppresses the Gibbs phenomenon 
ringing due to the binary edges of the world map. We then 
applied the filter ( [59] ) to obtain the result in Figure [5jc). The 
lowpass filtering of the 5-tap filter is evident. We then validated 
the theory further by using a directional filter. The butterfly 
filter is defined by McEwan et al (19] as follows: for (x^y) 
in the plane R^, let 



h{x, y) = xe 



(60) 



This filter may be mapped to the sphere using the co-latitude 
j3 and longitude a as 



h{P,a) = [tan(/3/2)cosa] 



(61) 



We generate an FIR filter from this prototype by using a set 
of N = 144 evenly-spaced samples on an angular grid in (3, 
and a. To determine an FIR filter we need a third Euler angle 
7, in order to set the rotation. We set 7 = —a to counter the 
rotation around z by a, and let Rk, for k = 1 to N, be the 
rotation defined by Euler angles ak, Pk^Jk = —OLk- Hence, the 
filter is specified on the sphere by the input-output equation 



N-l 



di^) = KPk,(^k)f{Rku). 



(62) 



In the frequency domain, the filter is specified at each fre- 
quency i by the transfer function 



AT-l 



H{1)= Y,KPk.o^k)D^{Ru). 



(63) 



/c=0 



We use the transfer function in our implementation, construct- 
ing H{t) for ^ = 0, . . . , 63. This filter is easily dilated by 
scaling each angle with a constant A as follows: Pk ^ 
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(a) 



(b) 




Fig. 4: (a) Shows the impulse function on a sphere reconstructed from its 64-th order spherical harmonic coefficients; (b) shows 
the impulse response of the 5-pt lowpass filter, with a central lobe and four sidelobes clearly visible; (c) shows the magnitude 
of the frequency response for the same filter, normalized by the magnitude spectrum of a single impulse to better illustrate its 
lowpass characteristic. 






(a) (b) (c) 

Fig. 5: Part (a) shows a binary world map, part (b) shows the reconstruction by spherical harmonics with bandwidth 63, and 



part (c) shows the reconstruction filtered by the 5-point lowpass filter in eq. (58), with the smoothing effect clearly visible. 



Qf/e ^ Aa/e. This relative ease of dilation suggests multi- 
resolution implementations, though we do not have space to 
pursue that here. Dilation may also be accomplished through 
stereographic projection d 12112911 or through harmonic scal- 
ing |32|. Figure [6] shows the results of filtering the world map 
with the FIR approximation to the butterfly, and the 90-degree 
rotation of it using y in place of x in ( [6Q| . Also shown are 
the results using a dilated filter with A = 2. 

In our third experiment, we extended the FIR filtering 
method proposed in this paper beyond the sphere S'^ to a 
larger set of two-dimensional manifolds. Let / : 5^ ^ 
be mapping with components denoted fx, fy, and fz- The 
components define the x, y, z coordinates of points on the 
surface of the manifold. Each of the three coordinate functions 
may be expanded using spherical harmonics. For each £, let 
denote the 3 x (2^ + 1) -dimensional matrix whose 3 rows 
contain the spherical harmonic coefficients of the correspond- 
ing coordinate function. Then the manifold equivalent of ( pT| ) 
is 



(64) 



£=0 



This expansion is referred to in the computer vision literature 
as SPHARM, the capital letters signifying that it is a 3- 



D version of spherical harmonics fT|. FIR filtering may be 
performed on SPHARM coefficients with a simple extension 
of ( [47] ), resulting in the equation 



(65) 



Note that in the above equation, Q and are 3 x (2^ + 1), 
while the transfer function matrix H{t) is 2^ + 1 -dimensional. 
Consequently, every filter developed for spherical harmonics 
may also be applied to SPHARM using exactly the same 
transfer function. Figure |7] shows the result of filtering the 
expansion of a cortical surface using the lowpass filter of ([59]), 
as well as the result of filtering with the product H{£)H{i) 
that is equivalent to cascaded filtering. 

In our experiments, computation time is dominated by the 
IRF, which takes approximately 16 seconds on a desktop 
computer with quad core 2.6GHz CPU for a L = 64 degree 
expansion of the world map data. Other steps, such the filtering 
operation ( [47] ), are negligible in comparison. The IRF may be 
replaced by the fast spherical harmonic transform 1 20 1 in cases 
of data where the Gibbs phenomenon caused by discontinuities 
is not a concern. 
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(a) 



(b) 





(c) 



(d) 



Fig. 6: Filtered versions of the world map. Part (a) shows the result of filtering using the directional butterfly filter (60), 
which acts as an edge detector, and part (b) shows the result with the orthogonally-oriented filter, which detects edges in a 
perpendicular direction. Part (c) shows the sum of absolute values of filter outputs in (a) and (b), illustrating that parts (a) and 
(b) form complementary aspects of an edge detector, and part (d) shows the output when the filter is dilated by a factor of 
two, showing a blurrier output as expected. 




\9 




(a) 



(b) 



(c) 



Fig. 7: Filtered versions of cortical surface. The manifold in part (a) is filtered in the SPHARM domain using the lowpass 



filter in eq. (58), giving the result shown in (b), and with a cascade of the same filter applied twice, giving a twice smoothed 
result as shown in part (c). 
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VI. Summary and Conclusions 

One of the key developments in the research leading to 
this paper is relatively simple: writing the spherical harmonic 
expansion in vector form ( pTj ). Once the expansion is seen 
in terms of vectors, the definition of phase, and its properties 
in relation to magnitude become tractable. It also becomes 
clear that to obtain a new set of vectors in a linear fashion, 
it is necessarily to multiply the existing set with compatible 
matrices. The vector-matrix form of the filtering described in 



Theorem 4.1 shows how such filtering is constructed. We have 
shown the results of experiments that validate the theory using 
both spherical and manifold data. 

One conclusion we reach in this work is that phase is more 
important than magnitude for spherical harmonics, and that 
a study of phase leads to better understanding of filtering. A 
second conclusion is that it is important to use vector-matrix 
notation, reducing the clutter of indices, to understand the 
operations involved in filtering. Finally, we conclude that FIR 
filtering on the sphere may be studied using the same formal 
methods of magnitude-phase analysis and transfer function 
derivation that serve so well in the familiar real-line case. 
This last conclusion suggests many opportunities for future 
work, to translate FIR filtering theory, including filter design 
and filter banks, to the spherical case using methods similar 
to those employed in this paper. 

Appendix A 

We prove ( [40| . Substituting from ( pTj ) into ([39]), and using 
the rotation property ([14]), we obtain 

9{R) = Y.Y.I FlY,{uyHkDk{R)Yk{u)du. (66) 
i k -^^^ 

Since the integrand is a scalar, we apply the trace and use its 
linearity and cyclic permutation invariance to write 



^(^) = 1^1] Trace / Ye{uy Ak{R)Yk{u)duF; 



(67) 



Here Ak{R) = Hj.Dj.{R) is a row-vector. Due to orthonor- 
mality ([6]), the integral vanishes if ^ ^ k, and for £ = k 
becomes the column vector Aj . Consequently, we have 

g{R) = ^ Trace [D^iRy F;] (68) 

i 

The result now follows by noting that Trace (X) = 
Trace(X^), and comparing with eq. i34i. 

Appendix B 



We prove ([44]). From (32) we have 



G{£) 



50(3) 



g{V)D,{V)UV. 



Substituting from (43), we get (all integrals are over S0(3)) 
that 



G{i) 



h{R)f{R-^V)dR 



Di{V)UV 



Reversing the order of integrals, we obtain 

G{i) = J h{R) J f{R-^V)Di{V)UV dR. 

Let U = R~^V, so that RU = V. Note that dV is a Haar 
measure on a compact group, so the integral is invariant to 
left shift (see for example ifTTi ). Consequently. 



G{1) = / h{R) 



f{U)D^{RU)UU 



dR. 



Now Di{RUy = Di{UyDi{R)\ so that we obtain 



G(^) = / h{R) 



Di{R)UR. 



f{U)D,{U)UU 
The term in brackets is F{i) by ( [32] ), and hence ( |44] ) follows. 

G{i)= [ h{R)F{i)Di{R)UR = F{£)H{i). 
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